%%html
<style>
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunsx.otf');
font-weight: bold;
}
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunsi.otf');
font-style: italic, oblique;
}
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunbxo.otf');
font-weight: bold;
font-style: italic, oblique;
}
.text_cell { font-family: "Computer Modern", "Palatino", "Palatino Linotype", serif; }
.code_cell:first-child { display: none; }
h1, h2, h3, h4 { text-align: center; }
</style>
%%capture
# Set-up and formatting
%pylab inline
from scipy.integrate import odeint
import scipy.signal as signal
# picture configuration
matplotlib.rcParams["savefig.dpi"] = 300
matplotlib.rcParams["savefig.bbox"] = "tight"
from IPython.core.interactiveshell import InteractiveShell
%config InlineBackend.figure_format = "retina" # higher-def plots
We begin by defining the constants for our system. The equations of motion for the triple pendulum are \begin{equation} \begin{pmatrix} \ddot\phi_1 \ \ddot\phi_2 \ \ddot\phi_3 \end{pmatrix} = \begin{pmatrix} 3 & 2\cos(\phi_1-\phi_2) & \cos(\phi_1-\phi_3) \ 2\cos(\phi_1-\phi_2) & 2 & \cos(\phi_2-\phi_3) \ \cos(\phi_1-\phi_3) & \cos(\phi_2-\phi_3) & 1 \ \end{pmatrix}^{-1} \begin{pmatrix} -\dot\phi_3^2\sin(\phi_1-\phi_3) -2\dot\phi_2^2\sin(\phi_1-\phi_2) -\frac{3g}{l}\sin\phi_1 \ -\dot\phi_3^2\sin(\phi_2-\phi_3) + 2\dot\phi_1^2\sin(\phi_1-\phi_2) -\frac{2g}{l}\sin\phi_2 \ \dot\phi_2^2\sin(\phi_2-\phi_3) +\dot\phi_1^2\sin(\phi_1-\phi_3) -\frac{g}{l}\sin\phi_3 \end{pmatrix}, \end{equation} which as we can see depend only on the angles and angular velocities of the masses, as well as $g$ and $l$.
# Constants
g = 9.8 # m/s^2
l = 0.5 # m
# Mechanics of triple pendulum system.
# Passed into odeint function.
def triple_pendulum(r, t, g, l):
phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = r # unpack variables
# precalculate trig functions
cos_1_2 = cos(phi_1 - phi_2)
cos_1_3 = cos(phi_1 - phi_3)
cos_2_3 = cos(phi_2 - phi_3)
sin_1_2 = sin(phi_1 - phi_2)
sin_1_3 = sin(phi_1 - phi_3)
sin_2_3 = sin(phi_2 - phi_3)
# define arrays and vectors
A = array([[3, 2*cos_1_2, cos_1_3],
[2*cos_1_2, 2, cos_2_3],
[cos_1_3, cos_2_3, 1 ]])
B = array([-dphi_3**2 * sin_1_3 - 2*dphi_2**2 * sin_1_2 - 3*g*sin(phi_1)/l,
-dphi_3**2 * sin_2_3 + 2*dphi_1**2 * sin_1_2 - 2*g*sin(phi_2)/l,
dphi_2**2 * sin_2_3 + dphi_1**2 * sin_1_3 - g*sin(phi_3)/l])
ddphi = matmul(inv(A), B) # angular acceleration
return concatenate(([dphi_1, dphi_2, dphi_3], ddphi)) # tack on angular velocity
def solve_triple_pendulum(r0, tmax, l=l, dt = 0.001):
t = arange(0, tmax, dt)
solution = odeint(triple_pendulum, r0, t, args = (g, l))
phi_1 = solution[:, 0]
phi_2 = solution[:, 1]
phi_3 = solution[:, 2]
dphi_1 = solution[:, 3]
dphi_2 = solution[:, 4]
dphi_3 = solution[:, 5]
return t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3
We also create some helper functions to plot solutions to the triple pendulum system. In addition to viewing the three angles over time, we'd like to be able to see the trajectories the masses trace out, as well as a phase space plot for each mass.
# Plots phi_1, phi_2, and phi_3 over time.
def time_plot(solution, subtitle = ""):
t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = solution
deg_1 = phi_1 * 180 / pi
deg_2 = phi_2 * 180 / pi
deg_3 = phi_3 * 180 / pi
figsize(15, 7)
fig = figure()
plot(t, deg_1, "-g", label="First pendulum")
plot(t, deg_2, "--m", label="Second pendulum")
plot(t, deg_3, ":k", label="Third pendulum")
# titles and labels
title("Pendulum Angles over Time\n{}".format(subtitle))
xlabel("Time (s)")
ylabel("Angle (deg)")
# legend
legend()
return fig
# Plots the pendulum trajectory in x-y space. NOT scaled appropriately to length l
def trajectory_plot(solution, subtitle = ""):
t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = solution
x1 = sin(phi_1)
x2 = x1 + sin(phi_2)
x3 = x2 + sin(phi_3)
y1 = -cos(phi_1)
y2 = y1 - cos(phi_2)
y3 = y2 - cos(phi_3)
figsize(15, 7)
fig = figure()
axvline(x = 0, color = "#bbbbbb")#vertical light gray line
axhline(y = 0, color = "#bbbbbb")#horizontal light gray line
plot(x1, y1, "-g", label = "First Pendulum")
plot(x2, y2, "--m", label = "Second Pendulum")
plot(x3, y3, ":k", label = "Third Pendulum")
title("Pendulum Trajectories over Time\n{}".format(subtitle))
xlabel("x / l")
ylabel("y / l")
legend()
margins(y = 0.1)
return fig
# Plots the phase space trajectories of the masses
def phase_space_plot(solution, subtitle = ""):
t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = solution
fig = figure()
title("Phase Space Trajectories of Masses\n{}\n\n".format(subtitle))
axis("off")
plot1 = fig.add_subplot(131)
axvline(x = 0, color = "#bbbbbb")# vertical light gray line
axhline(y = 0, color = "#bbbbbb")# horizontal light gray line
plot(phi_1, dphi_1, "k")
title("First Pendulum")
xlabel("Angle")
ylabel("Angular Velocity (1/s)")
plot2 = fig.add_subplot(132)
axvline(x = 0, color = "#bbbbbb")# vertical light gray line
axhline(y = 0, color = "#bbbbbb")# horizontal light gray line
plot(phi_2, dphi_2, "k")
title("Second Pendulum")
xlabel("Angle")
plot3 = fig.add_subplot(133)
axvline(x = 0, color = "#bbbbbb")# vertical light gray line
axhline(y = 0, color = "#bbbbbb")# horizontal light gray line
plot(phi_3, dphi_3, "k")
title("Third Pendulum")
xlabel("Angle")
tight_layout()
return fig
We want to check that the triple_pendulum function accurately captures the dynamic of the system. Since there is no analytic solution to the triple pendulum, we can do this by checking a few limiting cases, and checking that energy is conserved generally.
We first conduct a sanity check by ensuring that when $\phi_i(0)=0$ and $\dot\phi_i(0)=0$, the pendulum system remains motionless.
solution = solve_triple_pendulum([0, 0, 0, 0, 0, 0], tmax=10)
time_plot(solution, "No Initial Motion").savefig("figures/no_motion_time_sol.png")
We can also check that energy is conserved. We will start the bottom pendulum with a moderate amount of initial angular velocity ($\dot\phi_3(0)=7\ s^{-1}$), but leave the other intial conditions unchanged.
solution = solve_triple_pendulum([0, 0, 0, 0, 0, 7], tmax = 60)
time_plot(solution, "Moderate Intial Velocity of Third Pendulum").savefig("figures/moderate_velocity_time_sol.png")
trajectory_plot(solution, "Moderate Intial Velocity of Third Pendulum").savefig("figures/moderate_velocity_trajectory.png")
phase_space_plot(solution, "Moderate Intial Velocity of Third Pendulum").savefig("figures/moderate_velocity_phase_space.png")
# Energy conservation (cancelling common factor of m)
t, phi_1, phi_2, phi_3, dphi_1, dphi_2, dphi_3 = solution
U = -g*l*(3*cos(phi_1) + 2*cos(phi_2) + cos(phi_3))
T = (1/2)*l**2 * (3*dphi_1**2 + 2*dphi_2**2 + dphi_3**2 + 4*dphi_1*dphi_2*cos(phi_1-phi_2)
+ 2*dphi_2*dphi_3*cos(phi_2-phi_3) + 2*dphi_1*dphi_3*cos(phi_1-phi_3))
E = T + U
efig = figure()
plot(t, T, ":", label="Kinetic Energy")
plot(t, U, "--", label="Potential Energy")
plot(t, E, "-", label="Total Energy")
title("Energy Conservation")
xlabel("Time (s)")
ylabel("Energy/mass (J/kg)")
legend()
efig.savefig("figures/energy_conservation.png");
We see that while kinetic and potential energy vary over time, total energy remains constant.
We would like to be able to know more about a given solution than we can learn from studying its phase space plots—is the motion periodic? If so, what is the period? Are there multiple components of the system's frequency? Is there a natural frequency for the system?
We can use the Fourier transform, which transforms the angle-over-time signal into frequency space, to help us answer these questions. The numpy package comes with Fast Fourier Transform (FFT) functions to perform the discrete Fourier transform on a given array. In order to achieve sufficient frequency resolution, we must first simulate the system over a larger time interval.
When $\dot\phi_3=7$ initially, the system has a periodic solution. We can confirm this by looking at the Fourier transform of each angle.
def frequency_plot(solution, max_freq = 2):
t, phi_1, phi_2, phi_3, _, _, _ = solution
# take the Fourier transforms
fft_1 = fft.rfft(phi_1)
fft_2 = fft.rfft(phi_2)
fft_3 = fft.rfft(phi_3)
# calculate our time step and use it to scale our frequency axis
dt = t[1] - t[0]
freqs = fft.rfftfreq(len(phi_1), dt)
# fft returns complex numbers, and so we take the absolute value of each
figsize(15, 7)
fig = figure()
plot(freqs, absolute(fft_1), "-g", label = "First Pendulum")
plot(freqs, absolute(fft_2), "--m", label = "Second Pendulum")
plot(freqs, absolute(fft_3), ":k", label = "Third Pendulum")
title("Fourier Transform of Pendulum Angles")
xlabel("Frequency (1/s)")
ylabel("Angle")
xlim(0, max_freq) # we don't care about high frequencies
legend()
return fig
solution = solve_triple_pendulum([0, 0, 0, 0, 0, 7], tmax = 1000, dt = 0.02)
frequency_plot(solution).savefig("figures/fft_periodic.png");
We see that this particular solution has two main frequency components, which tells us the motion of the system is periodic. We can find the frequency that corresponds to the main peak using the argmax function, and then convert that to a period.
# recalculate FFT
t, phi_1, phi_2, phi_3, _, _, _ = solution
# take the Fourier transforms
fft_1 = fft.rfft(phi_1)
fft_2 = fft.rfft(phi_2)
fft_3 = fft.rfft(phi_3)
# calculate our time step and use it to scale our frequency axis
dt = t[1] - t[0]
freqs = fft.rfftfreq(len(phi_1), dt)
# Find the frequency at which the FFT is at a maximum
freq_1 = freqs[argmax(fft_1)]
freq_2 = freqs[argmax(fft_2)]
freq_3 = freqs[argmax(fft_3)]
# Confirm all three pendula have the same frequency
print("Main frequencies are {}equal.".format("" if freq_1 == freq_2 == freq_3 else "not "))
# Calculate period
period = 2*pi / freq_1
print("Period of main frequency: {:.4f} s".format(period))
When $\dot\phi_3(0)=12$, the solution is chaotic. We see this as an absence of any peaks in the frequency spectrum except for a massive singularity at 0, which suggests aperiodic motion.
solution = solve_triple_pendulum([0, 0, 0, 0, 0, 12], tmax = 1000, dt = 0.02)
frequency_plot(solution).savefig("figures/fft_chaotic.png");
We can explore the conditions under which the triple pendulum system becomes chaotic by making a contour plot of the FFT for each initial condition. Unlike the damped, driven pendulum, there is no natural frequency at which to sample the solution values, so we cannot make a traditional bifurcation diagram. Instead, we will plot the frequency spectrum for a range of initial angular velocities ($\dot\phi_3$), using color to denote the value of the Fourier transform at each frequency.
# Generate data for contour plot
# params: the index (or indices) of the parameter in the r0 vector which to plot over.
# For example, to simulate changing the initial angle of the 2nd pendulum from 0 to pi radians, we would call
# contour_data(1, 0, pi, ...) since phi_2 is stored in index 1 of the r0 vector
# start: the starting value of the parameter
# end: the ending value of the parameter
# xSize: the number of intermediate paramter values. Will be the x-resolution
# ySize: the number of frequencies between 0 and max_freq. Will be the y-resolution
# fft_param: which parameter to perform the FFT on and to plot. Defaults to phi_3
# max_freq: Maximum frequency to save. Defaults to 2 1/s
# file: name of the file to save the data to. All data will be saved to the data/ folder
def contour_data(params, start, end, xSize, ySize, fft_param = 3, max_freq = 2, file="contour"):
initial_values = linspace(start, end, xSize)
t_max = ySize / max_freq
data = zeros((xSize, ySize))
for step, p_0 in enumerate(initial_values):
# set up initial conditions vector
r0 = zeros(6)
r0[params] = p_0
# run simulation
solution = solve_triple_pendulum(r0, tmax=t_max, dt=0.05)
fft_data = absolute(fft.rfft(solution[fft_param]))[ : ySize] # get FFT and discard high frequencies
# normalize to this solution's values, concatenate onto our results array
data[step, :] = fft_data / max(fft_data)
save("data/{}.npy".format(file), data)
def plot_contour(data, start, end, max_freq = 2):
# plot
figsize(12,7)
fig = figure()
imshow(transpose(data), extent=[start, end, 0, max_freq], origin="lower",
interpolation="spline36", aspect=0.6*(end-start)/max_freq, cmap="plasma") # 0=white, 1=black
ylabel("Frequency of Solution (1/s)")
return fig
# WARNING: will take a long time to run.
contour_data(5, 2, 20, 600, 250, fft_param=3, file="contour_3")
data_3 = load("data/contour_3.npy")
fig = plot_contour(log(data_3), 2, 20)
title("Solution Frequency versus Initial Angular Velocity of Third Pendulum")
xlabel("Initial Angular Velocity (rad/s)")
fig.savefig("figures/contour_3.png")
contour_data(4, 2, 20, 600, 250, fft_param=2, file="contour_2")
data_2 = load("data/contour_2.npy")
fig = plot_contour(log(data_2), 2, 20)
title("Solution Frequency versus Initial Angular Velocity of Second Pendulum")
xlabel("Initial Angular Velocity (rad/s)")
fig.savefig("figures/contour_2.png")
contour_data(3, 2, 20, 600, 250, fft_param=1, file="contour_1")
data_1 = load("data/contour_1.npy")
fig = plot_contour(log(data_1), 2, 20)
title("Solution Frequency versus Initial Angular Velocity of First Pendulum")
xlabel("Initial Angular Velocity (rad/s)")
fig.savefig("figures/contour_1.png")
contour_data([0,1,2], 0.0001, 3*pi/4, 60, 250, fft_param=3, file="contour_angles")
data_angles = load("data/contour_angles.npy")
fig = plot_contour(log(data_angles), 0, 3*pi/4)
title("Solution Frequency versus Initial Angle of all Pendula")
xlabel("Initial Angle (rad)")
fig.savefig("figures/contour_angles.png")
(continued in Part II)